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Abstract One main issue, when numerically integrating autonomous Hamiltonian 
systems, is the long-term conservation of some of its invariants, among which the 
Hamiltonian function itself. Recently, a new class of methods, named Hamiltonian 
Boundary Value Methods (HBVMs) has been introduced and analysed [4|, which are 
able to exactly preserve polynomial Hamiltonians of arbitrarily high degree. We here 
study a further property of such methods, namely that of having, when cast as a 
Runge-Kutta method, a matrix of the Butcher tableau with the same spectrum (apart 
from the zero eigenvalues) as that of the corresponding Gauss-Legendre method, in- 
dependently of the considered abscissae. Consequently, HBVMs are always perfectly 
A-stable methods. This, in turn, allows to elucidate the existing connections with clas- 
sical Runge-Kutta collocation methods. 
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1 Introduction 

Hamiltonian problems are of great interest in many fields of application, ranging from 
the macro-scale of celestial mechanics, to the micro-scale of molecular dynamics. 
They have been deeply studied, from the point of view of the mathematical analysis, 
since two centuries. Their numerical solution is a more recent field of investigation, 
which has led to define symplectic methods, i.e., the simplecticity of the discrete 
map, considering that, for the continuous flow, simplecticity implies the conservation 
of H(y). However, the conservation of the Hamiltonian and simplecticity of the flow 
cannot be satisfied at the same time unless the integrator produces the exact solution 
(see |9| page 379]). More recently, the conservation of energy has been approached 
by means of the concept of the discrete line integral, in a series of papers 11 llll2l[T3l 
[T4l[T5l . leading to the definition of Hamiltonian Boundary Value Methods (HBVMs) 
EEHUO, which is a class of methods able to preserve, for the discrete solution, poly- 
nomial Hamiltonians of arbitrarily high degree (and then, a practical conservation of 
any sufficiently differentiable Hamiltonian). In more details, in (4), HBVMs based on 
Lobatto nodes have been analysed, whereas in J5] HBVMs based on Gauss-Legendre 
abscissae have been considered. In the last reference, it has been actually shown that 
both formulae are essentially equivalent to each other, since the order and stability 
properties of the method turn out to be independent of the abscissae distribution, and 
both methods are equivalent, when the number of the so called silent stages tends to 
infinity. In this paper this conclusion if further supported, since we prove that HB- 
VMs, when cast as Runge-Kutta methods, are such that the corresponding matrix of 
the tableau has the nonzero eigenvalues coincident with those of the corresponding 
Gauss-Legendre formula (isospectral property of HBVMs). 

This property can be also used to further analyse the existing connections between 
HBVMs and Runge-Kutta collocation methods. 

With this premise, the structure of the paper is the following: in Section[2]the basic 
facts about HBVMs are recalled; in Section [3] we state the main result of this paper, 
concerning the isospectral property; in Section [4] such property is further general- 
ized to study the existing connections between HBVMs and Runge-Kutta collocation 
methods; finally, in Section[5]a few concluding remarks are given. 

2 Hamiltonian Boundary Value Methods 

The arguments in this section are worked out starting from the arguments used in J4] 
13 to introduce and analyse HBVMs. We consider canonical Hamiltonian problems 
in the form 



where J is a skew-symmetric constant matrix, and the Hamiltonian H (y) is assumed 
to be sufficiently differentiable. The key formula which HBVMs rely on, is the line 
integral and the related property of conservative vector fields: 



y = JVH{y), y(t )=y €R 2m , 



(2.1) 




(2.2) 
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for any y\ € Mr", where o is any smooth function such that 

o(t )=y Ql a(t +h)=y l . (2.3) 

Here we consider the case where a(t) is a polynomial of degree s, yielding an ap- 
proximation to the true solution y(t) in the time interval [to, to +h]. The numerical 
approximation for the subsequent time-step, y\, is then defined by ( 12.3b . 
After introducing a set of s distinct abscissae, 

0<c u ...,c s <l, (2.4) 

we set 

Yi = a(t + Cih), i=l,...,s, (2.5) 

so that a if) may be thought of as an interpolation polynomial, interpolating the fun- 
damental stages Yi, i= 1, . . . ,s, at the abscissae ( 12.4l i. We observe that, due to ( 12.31 ), 
a(t) also interpolates the initial condition yo. 

Remark 1 Sometimes, the interpolation at to is explicitly required. In such a case, 
the extra abscissa cq = is formally added to ( 12. 4[ . This is the case, for example, of 
a Lobatto distribution of the abscissae 

Let us consider the following expansions of <j(f) and a(t ) for t £ [to,to + h]: 

s s 

&{t + Th) = £ YjPj(z), ct(/ + tA) = y + h £ y, / Pj(x) dx, (2.6) 

where {Pj(t)} is a suitable basis of the vector space of polynomials of degree at most 
s — 1 and the (vector) coefficients {y 7 } are to be determined. We shall consider an 
orthonormal basis of polynomials on the interval [0, 10 i.e.: 

/ Pi(t)Pj(t)dt = fy, i.j 1 v. (2.7) 

Jo 

where is the Kronecker symbol, and P,(f) has degree i— 1. Such a basis can be 
readily obtained as 

P,(0 = \ / 2^T^-i(0 5 i = l,...,J, (2.8) 

with Pj-\(t) the shifted Legendre polynomial, of degree / — 1, on the interval [0,1]. 

Remark 2 From the properties of shifted Legendre polynomials (see, e.g., STj or the 
Appendix in 13$), one readily obtains that the polynomials {Pj(t)} satisfy the three- 
terms recurrence relation: 

Pi(t) = l, P 2 (t)=V3(2t-l), 

1 The use of an arbitrary polynomial basis is also permitted and has been considered in the past (see for 
example 1151121 ). however as was shown in 0, among all possible choices, the Legendre basis turns out 
to be the optimal one. 
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We shall also assume that H(y) is a polynomial, which implies that the integrand 
in ( 12.21 i is also a polynomial so that the line integral can be exactly computed by 
means of a suitable quadrature formula. It is easy to observe that in general, due to 
the high degree of the integrand function, such quadrature formula cannot be solely 
based upon the available abscissae {c,}: one needs to introduce an additional set of 
abscissae {d\, . . . ,c r }, distinct from the nodes {c,}, in order to make the quadrature 
formula exact: 

/ 1 cr(*o + Th) T VH(a(t + th))dT = (2.9) 
Jo 

s r 

£ Pi6(t + Ci h) T WH(a(t + ah)) + £ Pi6(to + Cih) T WH{a(t Q + c t h)), 

1=1 !=1 

where J3,-, i = 1 , . . . , s, and /3,-, i = 1 , . . . , r, denote the weights of the quadrature formula 
corresponding to the abscissae {c,} U {c,}, i.e., 




(2.10) 



Remark 3 In the case considered in the previous Remark\l\ i.e. when cq = is for- 
mally added to the abscissae ( 12.41 ). the first product in each formula in H2.10[ ranges 
from j = to s. Moreover, also the range of {j3,} becomes i = 0, 1, ... ,s. However, 
for sake of simplicity, we shall not consider this case further. 

Definition 1 The method defined by the polynomial <j(t), determined by substituting 
the quantities in (12.61 l into the right-hand side of (12.91 , and by choosing the unknown 
coefficient {yf\ in order that the resulting expression vanishes, is called Hamiltonian 
Boundary Value Method with k steps and degree s, in short HB"VM(k,s), where k = 
s + r Ml?- 

According to fl4l . the right-hand side of j2.9i is called discrete line integral as- 
sociated with the map defined by the HBVM(fc,s), while the vectors 

Yi = a(tQ + Cih), i = l,...,r, (2.11) 

are called silent stages: they just serve to increase, as much as one likes, the degree 
of precision of the quadrature formula, but they are not to be regarded as unknowns 
since, from ( 12.6b and (12. lit , they can be expressed in terms of linear combinations of 
the fundamental stages (12.5b . 

In the sequel, we shall see that HBVMs may be expressed through different, 
though equivalent, formulations: some of them can be directly implemented in a com- 
puter program, the others being of more theoretical interest. 
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Because of the equality ( 12.91 ). we can apply the procedure described in Defini- 
tion Q] directly to the original line integral appearing in the left-hand side. With this 
premise, by considering the first expansion in (12.6b . the conservation property reads 



£yj [ l Pj(T)VH(a(t + Th))dT = Q, 

7=1 J0 

which, as is easily checked, is certainly satisfied if we impose the following set of 
orthogonality conditions: 

Yj= [ Pj{x)JVH{a{t + xh))dx, j = l,...,s. (2.12) 



Then, from the second relation of (12.61 ) we obtain, by introducing the operator 

L(f;h)a(t + ch)= (2.13) 

°(!o) + h £ fPj{x)dx f Pj(T)f(a(t + vh))dx, c e [0, l], 
Jo Jo 

that o is the eigenfunction of L(JVH;h) relative to the eigenvalue X = 1: 

o=L(TVH;h)o. (2.14) 

Definition 2 Equation (12.741 1 is the Master Functional Equation defining a ^j. 

Remark 4 From the previous arguments, one readily obtains that the Master Func- 
tional Equation \2.14\ characterizes HBVM(k 1 s) methods, for all k>s. Indeed, such 
methods are uniquely defined by the polynomial O, of degree s, the number of steps k 
being only required to obtain an exact quadrature formula (see ( 12.91 )). 

To practically compute a, we set (see ( 12.51 ) and ( 12.61 )) 

s 

Y i = o(to + c i h)=y + hY i aijYj, i=l,...,s, (2.15) 

7=1 

where 

dij = J Pj (x)6x, ij = l,...,s. 

Inserting ( 12.121 ) into ( 12.15b yields the final formulae which define the HBVMs class 
based upon the orthonormal basis {Pj}: 

Y i =y Q + hf j a i j f 1 Pj(T)JVH(a(t + Th))dT, i = l,...,s. (2.16) 

7=1 J ° 

For sake of completeness, we report the nonlinear system associated with the 
HBVM(A:,5) method, in terms of the fundamental stages {T,} and the silent stages 
{Yj} (see d2.1 j} ), by using the notation 

f{y)=JVH{y). (2.17) 
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In this context, it represents the discrete counterpart of ( 12.161 ). and may be directly 
retrieved by evaluating, for example, the integrals in (12.16b by means of the (exact) 
quadrature formula introduced in ( 12.91 i: 



From the above discussion it is clear that, in the non-polynomial case, supposing 
to choose the abscissae {c,} so that the sums in ( 12.18b converge to an integral as 
r = k — s —> °°, the resulting formula is ( 12.16b . 

Definition 3 Formula ( 12. 761 ) is named °°-HBVM of degree s or HBVM(°°,s) 

This implies that HBVMs may be as well applied in the non-polynomial case 
since, in finite precision arithmetic, HBVMs are undistinguishable from their limit 
formulae ( 12.16b . when a sufficient number of silent stages is introduced. The aspect 
of having a practical exact integral, for k large enough, was already stressed in ||2][4] 



On the other hand, we emphasize that, in the non-polynomial case, ( 12.161 ) be- 
comes an operative method only after that a suitable strategy to approximate the 
integrals appearing in it is taken into account. In the present case, if one discretizes 
the Master Functional Equation d2.13b - d2.14b . HB VM(fc, s) are then obtained, essen- 
tially by extending the discrete problem d2.18l ) also to the silent stages (12. lib . In order 



to simplify the exposition, we shall use (12.17b and introduce the following notation: 

{x i } = {c i }U{6 i }, {»,} = {ft} U {ft}, 
yi = a{t + Tih), ft = f{a{t + tih)), i = l,...,k. 
The discrete problem defining the HBVM(A:,.s) then becomes, 




(2.18) 



mm. 




k. 



(2.19) 



By introducing the vectors 



y = (y T 1 ,...,y T k ) T , 



e 



(i.-.ifeK', 



and the matrices 



Q =diag(fi)i,...,%), 



•kxs 



(2.20) 



whose (/, j)th entry are given by 





(2.21) 



we can cast the set of equations d2.19b in vector form as 

y = e <g> y + h{j? s @>jQ. ) <g> W(y) , 
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with an obvious meaning of /(y). Consequently, the method can be regarded as a 
Runge-Kutta method with the following Butcher tableau: 

















COi . . . COk 



(2.22) 



Remark 5 We observe that, because of the use of an orthononnal basis, the role of 
the abscissae {c,} and of the silent abscissae {c,} is interchangeable, within the set 
{t ; }. This is due to the fact that all the matrices J" s , £P S , and Q depend on all the 
abscissae {t,}, and not on a subset of them, and they are invariant with respect to the 
choice of the fundamental abscissae {c ( -}. 

In particular, when a Gauss distribution of the abscissae {tj , . . . , T&} is considered, 
it can be proved that the resulting HBVM(£,s) method [5|: 

- has order 2s for all k > s\ 

- is symmetric and perfectly A-stable (i.e., its stability region coincides with the 
left-half complex plane, C~ JD); 

- reduces to the Gauss-Legendre method of order 2s, when k = s; 

- exactly preserves polynomial Hamiltonian functions of degree v, provided that 

vs 

k>—. (2.23) 
- 2 

Additional results and references on HBVMs can be found at the HBVMs Home- 
page Q. 



3 The Isospectral Property 



We are now going to prove a further additional result, related to the matrix appearing 
in the Butcher tableau (12.22b . corresponding to HBVM(£,s), i.e., the matrix 



tkxk 



k>s, 



(3.1) 



whose rank is s. Consequently it has a (k — s)-fold zero eigenvalue. In this section, 
we are going to discuss the location of the remaining s eigenvalues of that matrix. 

Before that, we state the following preliminary result, whose proof can be found 
in iflOl Theorem5.6 on page 83]. 

Lemma 1 The eigenvalues of the matrix 

'•• 



x s = 



\ 







(3.2) 
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with 



2y/{2j +1)(2;-1)' 



7>1, 



(3.3) 



coincide with those of the matrix in the Butcher tableau of the Gauss-Legendre method 
of order 2s. 



We also need the following preliminary result, whose proof derives from the prop- 
erties of shifted-Legendre polynomials (see, e.g., Q] or the Appendix in J4)). 



Lemma 2 With reference to the matrices in \2.20\ - K2.21\ , one has 

where 



4i o 



x s = 



\ 



with the £j defined by ( li.il ). 



The following result then holds true. 

Theorem 1 (Isospectral Property of HBVMs) For all k> s and for any choice of 
the abscissae {t,} such that the quadrature defined by the weights {ft);} is exact for 
polynomials of degree 2s — 1, the nonzero eigenvalues of the matrix A in ( li.il ) coincide 
with those of the matrix of the Gauss-Legendre method of order 2s. 

Proof For k = s, the abscissae {t,} have to be the s Gauss-Legendre nodes, so that 
HB VM(j, s) reduces to the Gauss Legendre method of order 2s, as outlined at the end 
of Section |2] 

When k > s, from the orthonormality of the basis, see ( 12.7b . and considering that 
the quadrature with weights {»,} is exact for polynomials of degree (at least) 2s — 1, 
one easily obtains that (see ( I2.201 i- (I2.21I >) 

&JQ& S+1 = (I S 0), 

since, for all i = l,...,s, and j = l,...,s+l: 

(£?JQ& S+1 ) = £ a>ePi(tt)Pj(t e ) = f Pi(t)Pj(t)6t = fy. 
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By taking into account the result of Lemma[2l one then obtains: 

A^ s+l = J? s ^jQ<? s+l = J s (7 V 0) = & S+1 Z S (/,. 0) = &> s+l (X s 0) 



(\ -«1 

4i o 



V" 



(3.4) 



0/ 



with the defined according to ( 13.3b . Consequently, one obtains that the columns 
of & s +\ constitute a basis of an invariant (right) subspace of matrix A, so that the 
eigenvalues of X s are eigenvalues of A. In more detail, the eigenvalues of X s are 
those of X s (see ( 13.21 )) and the zero eigenvalue. Then, also in this case, the nonzero 
eigenvalues of A coincide with those of X s , i.e., with the eigenvalues of the matrix 
defining the Gauss-Legendre method of order 2s. □ 



4 HBVMs and Runge-Kutta collocation methods 



By using the previous results and notations, now we further elucidate the existing 
connections between HBVMs and Runge-Kutta collocation methods. We shall con- 
tinue to use an orthonormal basis {Pj}, along which the underlying extended collo- 
cation polynomial a if) is expanded, even though the arguments could be generalized 
to more general bases, as sketched below. On the other hand, the distribution of the 
internal abscissae can be arbitrary. 

Our starting point is a generic collocation method with k stages, defined by the 
tableau 



(4.1) 

















CO!... CO k 



where, for i, j — l,...,k: 

stf = (ay) = 



tj{x)Ax I , ®j = / £j(x)dx 



being the j'th Lagrange polynomial of degree k — 1 defined on the set of abscis- 
sae {t,}. 

Given a positive integer s < k, we can consider a basis {pi(x),...,p s (x)} of the 
vector space of polynomials of degree at most 5—1, and we set 



( Pl(tl) P2(Tl) ••• Ps{l\)\ 



\P\M P2(*k) ■■■ Ps{Tk)J 



(4.2) 



k x .v 
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(note that tP s is full rank since the nodes are distinct). The class of Runge-Kutta 
methods we are interested in is defined by the tableau 



(4.3) 



<0] 



(Ok 



where £2 = diag(oi , . . . , (0%) (see ( 12.201 )) and A s = diag(r/i , . . . , T7 S ); the coefficients 
rjj, j = I,... ,s, have to be selected by imposing suitable consistency conditions on 
the stages {y,} (see, e.g., Q). In particular, when the basis is orthonormal, as we shall 
assume hereafter, then matrix & s reduces to matrix £P S in (12.20b - d2.2U . A s = I s , and 
consequently (14.3b becomes 



(4.4) 



(0\ 



(Ok 



We note that the Butcher array A has rank which cannot exceed s, because it is 
defined by filtering stf by the rank s matrix Q,. 

The following result then holds true, which clarifies the existing connections be- 
tween classical Runge-Kutta collocation methods and HBVMs. 

Theorem 2 Provided that the quadrature formula defined by the weights {©,-} is 
exact for polynomials at least 2s — 1 (i.e., the Runge-Kutta method defined by the 
tableau \4. 41 ) satisfies the usual simplifying assumption B(2s)), then the tableau \4.4\ 
defines a HBVM(k,s) method based at the abscissae {t,}. 



Proof Let us expand the basis {P\ (t), . . 
j = l,...,k, defined over the nodes t,-, i = 1 , 



,P v (t)} along the Lagrange basis {£ 7 (t)}, 
...fc 



= 1, 



.,s. 



It follows that, for i = 1 , . . . , k and j = 1 , . . . , s: 



f'p j {x)dx= YPj{Xr) f' l r {x)dx- 
JQ ,fi JO 



that is (see d220l - d23B and KT\ ). 



r=l 



Ob, 



(4.5) 



By substituting ( 14.51 ) into ( 14.4b . one retrieves that tableau ( 12.22b . which defines the 
method HBVM(£,s). This completes the proof. □ 

The resulting Runge-Kutta method ( 14.4b is then energy conserving if applied to 
polynomial Hamiltonian systems ( 12.1b when the degree of H(y), is lower than or 
equal to a quantity, say v, depending on k and s. As an example, when a Gaussian 
distribution of the nodes {t, } is considered, one obtains ( 12.23b . 
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Remark 6 (About Simplecticity) The choice of the abscissae {t\ , . . . , x k } at the Gaus- 
sian points in [0, 1] has also another important consequence, since, in such a case, 
the collocation method H4. 1 i is the Gauss method of order 2k which, as is well known, 
is a symplectic method. The result of Theorem\2\then states that, for any s < k, the 
HBVM{k 1 s) method is related to the Gauss method of order 2k by the relation: 

where the filtering matrix (3^ S 3^J Q) essentially makes the Gauss method of order 
2k "work" in a suitable subspace. 

It seems like the price paid to achieve such conservation property consists in the 
lowering of the order of the new method with respect to the original one ( 14.11 ). Actu- 
ally this is not true, because a fair comparison would be to relate method ( I2.22I )-( I4.4I ) 
to a collocation method constructed on s rather than on k stages, since the resulting 
nonlinear system turns out to have dimension s, as shown in |4|. This computational 
aspect is fully elucidated in a companion paper (6), devoted to the efficient implemen- 
tation of HBVMs, where the Isospectral Property of the methods is fully exploited 
for this purpose. 



4. 1 An alternative proof for the order of HBVMs 

We conclude this section by observing that the order 2s of an HBWM(k,s) method, 
under the hypothesis that ( 14. U satisfies the usual simplifying assumption B(2s), i.e., 
the quadrature defined by the weights {©,} is exact for polynomials of degree at least 
2s— 1 , may be stated by using an alternative, though equivalent, procedure to that 
used in the proof of [4, Corollary 2] (see also J3] Theorem 2]). 

Let us then define the k x k matrix = &> k (see d2.20b - (l2.211 i) obtained by 
"enlarging" the matrix & s with k — s columns defined by the normalized shifted 
Legendre polynomials Pj{i), j = s + 1, . . . ,k, evaluated at {t,}, i.e., 



Pl(z k ) ... P k (T k ) 



By virtue of property B(2s) for the quadrature formula defined by the weights {ft),}, 
it satisfies 



OR 

This implies that & satisfies the property T(s,s) in ifTUl Definition 5. 10 on page 86], 
for the quadrature formula (ft),-, T/)* =1 . Therefore, for the matrix A appearing in ( 14.41 ) 
(i.e., ( 12.221 ). by virtue of Theorem|2]i, one obtains: 

&>- x A3? = &>- x $i3t ( Is _ ) = ( Xs \ , (4.6) 



O) \ o t 

where X s is the matrix defined in ( 13.41 ). Relation ( 14.6b and IfTUl Theorem 5. 1 1 on page 
86] prove that method ( 144b (i.e., HBVM(£,s)) satisfies C(s) and D(s- 1) and, hence, 
its order is 2s. 
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5 Conclusions 

In this paper, we have shown that the recently introduced class of HB VMs (k, s), when 
recast as Runge-Kutta method, have the matrix of the corresponding Butcher tableau 
with the same nonzero eigenvalues which, in turn, coincides with those of the matrix 
of the Butcher tableau of the Gauss method of order 2s, for all k > s such that B(2s) 
holds. 

Moreover, HBVM(k,s) defined at the Gaussian nodes {ti , . . . , TjJ on the interval 
[0, 1] are closely related to the Gauss method of order 2k which, as is well known, is 
a symplectic method. 

An alternative proof of the order of convergence of HBVMs is also provided. 
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